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Abstract 

In a number of geophysical or planetological settings, including Earth's inner core, a sili- 
cate mantle crystallizing from a magma ocean, or an ice shell surrounding a deep water ocean 
- a situation possibly encountered in a number of Jupiter and Saturn's icy satellites -, a con- 
verting crystalline layer is in contact with a layer of its melt. Allowing for melting/freezing 
at one or both of the boundaries of the solid layer is likely to affect the pattern of convection 
in the layer. We study here the onset of thermal convection in a viscous spherical shell with 
dynamically induced melting/freezing at either or both of its boundaries. It is shown that the 
behavior of each interface - permeable or impermeable - depends on the value of a dimen- 
sional number V (one for each boundary), which is the ratio of a melting/freezing timescale 
over a viscous relaxation timescale. A small value of V corresponds to permeable boundary 
conditions, while a large value of V corresponds to impermeable boundary conditions. The 
linear stability analysis predicts a significant effect of semi-permeable boundaries when the 
number V characterizing either of the boundary is small enough: allowing for melting/freezing 
at either of the boundary results in the emergence of larger scale convective modes. The ef- 
fect is particularly drastic when the outer boundary is permeable, since the degree 1 mode 
remains the most unstable even in the case of thin spherical shells. In the case of a spherical 
shell with permeable inner and outer boundaries, the most unstable mode consists in a global 
translation of the solid shell, with no deformation. In the limit of a full sphere with permeable 
outer boundary, this corresponds to the "convective translation" mode recently proposed for 
Earth's inner core. As an example of possible application, we discuss the case of thermal 
convection in Enceladus' ice shell assuming the presence of a global subsurface ocean, and 
found that melting/freezing could have an important effect on the pattern of convection in 
the ice shell. 



1 Introduction 



The seismologically observed hemispherical asymm etry of the inner core [Tanaka and Hamaguchil . 
1997 : Niu and Wen . 2001 : Irving and Deussl 2011 1 has recently been interpreted as resulting from 
a high- viscosity mode of thermal convection, consisting in a translation of the inner core with 
melti ng on one hemisphere and solidification on the other Monnereau et al. , 2O10l : Alboussiere et al 



2010l | . This "convective translation" regime can exist because the boundary between the inner 
core and the outer core is a phase change interface, which means that deforming the inner core 
boundary (ICB) by internal stresses can induce melting or freezing. Melting occurs when the 
ICB is displaced outward, and crystallization occurs when the ICB is displaced inward, at a rate 
which depends on the ability of outer core convection to supply or evacuate the latent heat of 
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phase change. Because there is no deformation, and therefore no viscous dissipation, associated 
with it, the translation mode is dominant whenever phase change at the inner core boundary 
proceeds at a fast enough rate. 

The situation where a convective crystalline shell is in contact with its melt is encountered in 
a number of other geophysical or planetological problems, including convection i n a silicate man- 



tle crystallizing from below from a magma ocean, or from a basal magma ocean Labrosse et al 



20071 : iTJlvrova et all l2012j | , or convection in an ice shell surrounding a dee p water ocean, a situ 



ation possibly encountered in several of Jupiter and Saturn's icy satellites jKivelson et all [2000J; 



Spohn and Schubertl . 120031 : iTylerl . 120081 ] . If one of the boundary is impermeable, the translation 
mode predicted for a full sphere obviously cannot exist, but we might anticipate that allowing 
for phase ch ange at the other boundary will modify the pattern of convection and favor larger 
scale modes [Monnereau and Dubuffetl . 120021 ] . 

We will study here the onset of thermal convection in a uniformly heated spherical shell with 
boundary con ditions allowing for dynam ically induced melting or freezing at either or both of 
the boundary [Deguen et all lsubmitted| . The problem set-up and the system of equations are 
described in section O with an emphasis on the formulation of the boundary conditions. The 
steady basic solution of the system of equations is found in section [3l In section 0J we perform a 
linear stability analysis of the set of equations described in section [21 which allows to determine 
the pattern of the first unstable mode as a function of the shell outer-to-inner radius ratio and of 
two non-dimensional numbers describing the resistance to phase change at each boundary. The 
results and possible applications are discussed in section [5j 



2 Problem definition 

We consider a viscous solid spherical shell of outer radius R and inner radius jR, in contact 
with melt layers either above or below, or both (see Figure [1]). Superscripts "+" or "-" will 
be used for quantities taken at the outer or inner boundary, respectively. The solid shell has 
constant density p s , the layers below and above have densities p~ and respectively, and we 
note Ap + = /?+ — p s and Ap~ = p~ — p s . To insure long term mechanical stability of the solid 
layer, we must have p~ > p s > p+, or Ap + > and Ap~ < 0. The inner and outer boundaries 
are phase change interfaces, and melting and freezing can therefore occur when the interface is 
displaced by internal stresses. This will be described (section 12. It) with a parametrization of the 
relationship between the freezing or melting rate and the dynamic topograp hy of the interface , 



which has been dev e loped for t he describing convection in Earth's inner core [Alboussiere et al 



2010 : iDeguen et all Isubmitted | . 



To be consistent with the assumption of constant density p s , the acceleration of gravity g in 
the spherical shell is assumed to vary linearly with radius r, g = —g're r , where g' = dg/dr = 
g + /R = C st , which is relevant to situations where the depth dependence of the density is too 
small to have a significant effect on the mean gravity profile. While this is not true in a number 
of situations of geophysical interest (like in Earth's mantle), we will make this assum ption for 
two reasons: (i) it is (mathematically) the simplest configuration Chandrasekharl . 1 196 11 ] . and (ii) 
the case of the inner core, for which g is essentially linear in r, corresponds to the limit 7 — > of 
the problem discussed here. Considering a more general form for g is likely to give qualitatively 
similar results. 

The spherical shell is heated volumetrically at a rate p s c ps S (with S in K/s). The rheology is 
assumed to be Newtonian and temperature and pressure independent, with a constant viscosity 
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Figure 1: A sketch of the problem considered here. 



r\. Thermal convection in the spherical shell is then described by the conservation equations for 
mass, momentum, and entropy, which take the form 

Vu = 0, (1) 
= - Vp - a Ps G g + r/V 2 u, (2) 

<9G 

— + u • VG = kV 2 G + S, (3) 
under the Boussinesq approximation. 
2.1 Boundary conditions 

The rate of melting/freezing at each interface depends on the ability of convective motion in the 
melt layer to transport the heat absorbed or released by the phase change. Given a topography 
h(9, <j)) of the boundarjfl, the rate of erosion of the topography by melting or freezing is set by 
a balance between the rate of latent heat release or absorption, p s L dh/dt, with the convective 
heat flux on the melt side, which should scale as pic p iu'SQ, where L is the latent heat of melting, 
c p i the specific heat capacity of the melt, v! a typical velocity scale for convective motion in the 
melt layer, and 5®(9,4>) the difference of potential temperature between the boundary and the 
adjacent melt. The b oundary is assumed to rem ain very close to thermodynamic equilibrium 
(more justifications in Deguen et al. submitted! ) . and is therefore at the melting temperature 



Mefined here in reference to the isopotential surface which coincides on average with the boundary. 
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T m . The potential temperature variation 50 along the boundary results from the combined effect 
of the pressure dependency of T m and of the adiabat in the melt layer, so that a topography h 
induces a difference of potential temperature between the boundary and the melt layer given by 



50 



-(m p - m ad )pfg ± h 



(4) 



where m p = dT s /dP is the Clapeyron slope, and m a & = dT^/dP the adiabatic gradient in the 
melt layer. With this expression for 50, the heat balance described above gives 



U r — 



dh 

at 



(5) 



where the timescale for phase change, , is 
PsL 



pfc p i \m p - m ad | g ± u r 



(6) 



With m p = T m A/9 ± /(/9 s/ o^ = L) from the Clapeyron relation, Equation ([6]) can be rewritten as 



n 2 T 2 
Ps L 



Pi {Ap^CpiTrn (1 - m ad /m p )g ± u' 



(7) 



Assuming that the phase-change timescale and the viscous relaxation timescale t v = rj/ (| A/9 ± |g ,± i?) 
are both small compared to the dynamical timescale of the shell (overturn time), we can neglect 
dh/dt in Eq. ([5]), which gives the boundary condition 



h_ 



(8) 



The mechanical boundary conditions are tangential stress-free conditions and continuity of 
the normal stress at both boundary. Under the assumption of small topography, the stress-free 
tangential condition writes 



Tr9 = V 



d / ue \ ldu 
r ir — +- 



dr \ r 



Tr<f> = V 



d_ /«£ 

dr V r 



+ 



r 80 
1 du r 



r sin! 



(9) 
(10) 



at r = 7 and 1, where t t q and T r( j, are the (r, ff) and (r,4>) components of the deviatoric stress 
tensor r. Continuity of the normal stress at each boundary is written as 



■pI 



8u r 



0, 



(11) 



il h 



where [. . .] denotes the difference of a quantity across the boundary. When expanded around 
the mean position of the boundary, Eq. pip gives 



. i , du r 
-Ap ± g h - 2n—- + p=0 
or 



(12) 
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under the assumption that pressure fluctu ations on th e melt side are negligible compared to 
pressure fluctuations on the solid side [e.g. iRibd . 120071 ] . With h related to u r by Eq. (jSJ), Eq. 
(|12|) gives a boundary condition for u r only: 



2r?^+p = 0. 



(13) 



The topography h is an implicit variable of the problem, and can be calculated a posteriori from 
the radial velocity at the boundary. 



2.2 Non-dimensional set of equations 

The governing equations and boundary conditions are now made dimensionless using the thermal 
diffusion timescale k/R 2 , the outer radius R, k/R, ijk/R 2 and SR 2 /(Qk) as scales for time, 
length, velocity, pressure and potential temperature, respectively. Using the same symbols for 
dimensionless quantities, the system of equations (p3i3]) is then written as 



Vu 




0, 



-Vp + RaQr + V 2 u, 



~dt 



+ u-VG = V 2 6 + 6, 



where the Rayleigh number is defined as 
a Ps g+SR 5 



Ra 



6r/K 2 



(14) 
(15) 

(16) 



(17) 



The Rayleigh number defined here is based on the outer radius R, not the shell thickness (l — y)R. 
Also, note that the Rayleigh number used here is half that defined by IChandrasekharl 19611 ] . The 
dimensionless boundary conditions at r = 7 or 1 can be written 



G( 7 ) = e 7 , 9(1) = 0, 



d t uq \ 1 du r 
dr v r / r d9 



1 du r 



dr \ r 

±V ± u r + 2 



0, 



r sin V 

du r 
dr 



P = 0, 



(18) 
(19) 

(20) 



where the "phase change numbers" V + and V are defined as [Deguenl . 120121 : Dcgucn ct al., 
submitted ] 



(21) 



where is the timescale for erosion of a topography by melting or freezing, as defined in Eq. 
([6]), and t v = rj/(\ /S.p ± \g ± R) is the viscous relaxation timescale at the lengthscale R. The phase 
change numbers ~P + and V~ are measures of the resistance to phase change on each boundary. 
In the limit of infinite , the boundary condition (|20p reduces to the condition u r = 0, which 
corresponds to impermeable conditions. In contrast, when 7 ?± —¥ 0, Eq. (120p implies that the 
normal stress tends toward at the boundary, which corresponds to fully permeable boundary 
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conditions Monnereau and Dubuffet . 2002( |. The general case of finite gives boundary con- 



ditions for which the rate of phase change at the boundary (equal to u r ) is proportional to the 
normal stress induced by convection within the spherical shell. Note that we have defined V~ 
and V + using the absolute value of Ap 1 * 1 , so that both V and V + are positive. Because Ap 
is negative, this introduces a minus sign before V~ in the boundary condition (|20p for the inner 
boundary. 



With the assumptions made so far, the velocity field is known to be purely poloidal Ribc, 



20071 ]. and we introduce the poloidal scalar P defined such that u = V x V x (Pr). Taking the 



curl of the momentum equation (115j) gives 

RaL 2 e= (V 2 ) 2 L 2 P, (22) 
where the angular momentum operator L 2 is 

r2 1 9 ( ■ a 9 \ 1 92 d*\ 



Horizontal integration of the momentum equation Q15|) [Ribd . 120071 ] shows that, on both boundary. 



-p+| : (rV 9 P)=C7-. (24) 

Using this expression to eliminate p in the boundary condition (|20|) . and noting that u r = ^L 2 P, 
continuity of the normal stress at each boundary (equation ([20]) ) gives the following boundary 
condition for the poloidal scalar at r = 1 or 7: 

9 ( rV 2 P - -L 2 p) - ±V±— = C s \ (25) 



dr \ r 
while the stress-free conditions (119j) give 

°~ P + (L 2 - 2) ^ = C st , r = 7 or 1. (26) 



3 Steady basic solution 

The governing equations and boundary conditions presented in section I2T21 admit a steady solution 
(denoted by an overbar .7.) in which the velocity field is u = and the potential temperature 
field 6 is given by the steady state, conductive version of Eq. ([3]), which writes 

= V 2 9 + 6. (27) 

With 6>(r = 1) = 0, the general solution of Eq. (|27p is of the form 

^ 1 — a 9 . 

9 = a + r 2 , (28 

r 

where the constant a depends on the thermal boundary condition (imposed temperature or flux) 
at r = 7. The stability analysis could be carried out for the general potential temperature profile 
given by Eq. (I28p . but we will here consider only the case a = 1. This is mathematically simpler, 
and, in addition, will allow us to extrapolate easily the results to the case of Eart h's inner core , 



for which t he basic diffusive potential temperature profile is given by O = 1 — r [De guen et al 
submitted] . The potential temperature at r = 7 is 6(7) = 1 — 7 2 . 
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4 Linear stability analysis 



We now investigate the stability of the basic conductive state against infinitesimal perturbations 
of the temperature a nd velocity fields. The present analysis follows the analysis presented in 



Chandrasekharl 19611 ] (Chapter VI-60), where the stability analysis is treated in the case of im- 



permeable boundaries, which corresponds to the limit of infinite V~ and V + . The case of thermal 
convection in a full sphere with boundary conditions as described above, which corresponds to 
the limit 7 — >• of the problem considered here, has been treated in lDeguen et al. submitted| . 

The temperature field is written as the sum of the conductive temperature profile given by 
Eq. (|28j) and infinitesimal disturbances 0, Q(r,9,4>,t) = 6(r) + Q(r,9,4>,t). The velocity field 
perturbation is denoted by u(r, 9, (ft, t), and has an associated poloidal scalar P(r,0,(f>,t). We 
expand the temperature and poloidal disturbances in spherical harmonics, 

00 1 



@ = E E tT(r)Y l m (0, ( f>)e^, (29) 

1=0 m=-l 
oo I 

P = T, E iW^ m (M)e^, (30) 

1=1 m=—l 

where 07 is the growth rate of the degree I perturbations (note that since m does not appear in 
the system of equations, the growth rate is function of I only, not m). 

The only non-linear term in the system of equations is the advection of heat u- V0 in Equation 
(I16|) . which is linearized as 

u r — = -2ru r = -2L 2 P. (31) 
or 

The resulting linearized transport equation for the potential temperature disturbance is 

!L - V 2 J 8 = 2L 2 P + 6. (32) 
at J 

Using the decompositions (j29|) and (|30p . the linearized system of equations is then, for I > 1, 

Ratr = Vfpr, (33) 
(a l -V l )tr = 2l(l + l)pr, (34) 

with the stress-free boundary condition written as 

-JL + [I(Z + 1) - 2] 3L = 0, I>1, (35) 

with r = 1 or 7 on the upper or lower boundaries, and the boundary conditions derived from the 
continuity of the normal stress given by 

- [rV lP T - 21(1 + 1)?LJ = +1(1 + l)V+(t) P J-, (36) 

J; (rV W T - 21(1 + l)?pj = -1(1 + l)V-(t) P -^, (37) 
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at the outer and inner boundaries, respectively (note the different signs of the right-hand-side 
terms). The operator T>i is defined as 



v = d 2 | 2 d 1(1 + 1) 
dr 2 r dr r 2 

We expand the potential temperature perturbation t™(r) as 

j 

where the functions C/;(ayr) are defined as 

Cu(ai jr ) = J_^ + i)(ay 7 )J z+ i(a y r) - J l+ i (e*y7) J_( z+ i) (a«j r ) 



(38) 



(39) 



(40) 



Chandrasekharl . Il96ll ] . Here J/% denotes the Bessel function of the first kind of degree k, and the 



co nstants a/,- are the jth zeros of the function Cu(r). By construction, Cu(aij^f) = 0. As discussed 
by IChandrasekharl 19611 ] . the functions Cu(otjr) form an integral set of functions satisfying the 
orthogonality relation 



C u (aijr)Cu(ai k r)rdr = N l+lj 5 jk , 



(41) 



where 



Ni +y - M 



Writing the poloidal scalar perturbations pf 1 as 
p7(r) = Y, A l3PiA r ) 



(42) 



(43) 



and injecting the expansions of £p(r) and p™(r) given by Eq. (]39p and fj43f) in the momentum 
equation ([33]) , the functions pij are solutions of the equation 



Noting that 

Cuiaijr) 



V, 



-a 



equation 



Plj 



r J ^Jr 

has a general solution of the form 



RaCu(a^r)_ + y y+2 ^ r _ (I+a) + B i r -(J-i). 
at \fr v 1 6 4 



(44) 



(45) 



(46) 



The coefficients B\ 4 are determined by the boundary conditions at the inner and outer bound- 
aries of the shell, as explained in Appendix [Al 
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Injecting the above solution for py and the potential temperature expansion (f3"9"j) in the 



linearized heat equation (IM 



in Aij (see IChandrasekhar 



, we obtain after some manipulation an infinite set of linear equations 



19611 ] ). which admits a non trivial solutio n only if its determina nt is 



equal to zero. With our choice of basic state and g oc r, and following I Chandrasekharl [19611 ]. we 
find a characteristic equation of the form 



N, 



Ra 



Qkj 



o, 



where 



Qkj 



1(1 + 1)2 cr k 

denotes the determinant, and where the functions Q^j are defined as 

(i-l) 



(47) 



Cu(a k r) B{r l + B j 2 r l+2 + B ] 3 r~ {l+1) + B{i 



r^dr. 



(48) 



Solving Eq. (|4"7|) with the B\ A determined by the boundary conditions (Appendix [A]) gives the 
critical Rayleigh number Ra l c for a perturbation of degree I. The pattern of the first unstable 
modes can be calculated by solving the system in A\j for given V" , V + and Ra, which then 
allows to calculate the poloidal scalar pf 1 from equations (f43l) and (l46j) . 



5 Results and applications 

Figure [2] shows the critical Rayleigh number corresponding to the degree one mode as a function 
of V~ and V + for three configurations : (i) the outer boundary is impermeable (V + 3> 1) and 
V~ is varied from permeable to impermeable conditions; (ii) the inner boundary is impermeable 
(V~ 3> 1) and V + is varied from permeable to impermeable conditions; and (iii) V + = V~ , with 
boundary conditions varied from permeable to impermeable. For all three configurations, there 
is a marked change in the critical Rayleigh number at some transitional value of V~ or V + ', with 
the critical Rayleigh number being significantly smaller when V~ or V + are smaller than this 
transitional value, corresponding to permeable conditions. 

In what follow, we will focus on end-members cases, for which each boundary is either perme- 
able (V^ "C 1) or impermeable (V^ S> 1), which gives four end-member configuration. Figure [3] 
shows the critical Rayleigh number as a function of I and 7 for the four end-member cases. The 
pattern of the first unstable mode (as well as the second for the V~ <C 1, V + S> 1 cases) are 
shown in Figure [H The degree of the first unstable mode is shown in Figure [5] as a function of 7 
for the four end-member cases. Each end-member case is described below. 



5.1 (V + ,V ) S> 1 — impermeable boundaries 

Letting V + and V~ tend toward infinity, the problem tends toward the case of Rayl eigh-Benard 
conve ction in a spherical shell with impermeable stress free boundaries, as discussed in I Chandr asekhar 
196 1| . and will be used a s a reference case for the present study. The results found here are 



identical to that found by IChandrasekharl 1961 1 ( except that, a s explained above, the Rayleigh 
numbers shown here are half that found by I Chandrasekharl 19611 ] because of different definitions). 
The degree one mode is the first unstable mode for 7 smaller than ~ 0.23. The degree of the 
first unstable mode then increases rapidly when 7 is increased (Figure [5]). The corresponding 
wavelength is commensurate with the shell thickness 1 — 7: assuming a relationship of the form 
l c = a/(l — 7) + b (which, given that A c ~ l/l c when l c S> 1, is equivalent to A c ~ 1 — 7), least 
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square inversion of £ c (l) gives l c = 2.17/(1 — 7) — 1.35, which is shown as a grey dash-dotted 
line in Figure [5j The fit is indeed good, consistent with the assumption of a critical wavelength 
proportional to the layer thickness. 



5.2 (V + ,V ) < 1 - permeable inner and outer boundaries 

On the other extreme, when both boundaries are fully permeable, the first unstable mode is 
always the degree one mode (Figure [3] and [5]) , which takes the form of a solid translation of the 
spherical shell (Figure H] and Appendix [B]). The limit of a full sphere (7 = 0) corresponds to 



the "convective translat i on" m ode recently put forward for Earth's inner core Monnereau et al 



20ld : lAlboussiere et all l20ld | 



Since this mode consists of a pure translation, there is no deformation, and therefore no 
viscous dissipation in the inner core. This of course does not mean that this is a non-dissipative 
mode. There is viscous (and magnetic in the case of Earth's inner core) dissipation in the melt 
layer associated with the redistribution of the latent heat of phase change. The melt layer must 
provide mechanical work to account for the dissipation associated with the redistribution of the 
latent heat, which means that this mode of convection is ultimately limited by the vigor of 
convective motions in the melt layer. 

It can be shown (Appendix [Bj that the emergence of the translation mode requires that the 
quantity 



R v 



Ra 



V+ + 7 2 V- 



(49) 



is higher than a critical value which is a function of 7 only. The quantity V + +7 2 V~ is, save for a 
factor 1 + 7 2 , the boundary area weighted mean of V~ and V + . Figure [6] shows the critical value 
of R-p for the translational instability as a function of 7 calculated using Eq. (J72J) of Appendix 
El When 7 — > 0, the criti cal value tends toward the value of (Ra/T ,+ ) c = 175/2 = 87.5 found by 



Deguen et al.l submitted! for a full sphere. R-p then increases with 7 



T he limit 7 — > is relevant for Earth's inner core dynamics [Monnereau et al.l . l2010l : IAlboussiere et al 



20101 : iMizzon and Monnereau! . 120131 : iDeguen et all lsubmitted| . The case of a spherical shell with 



phase change at both boundaries might be relevant for the early dynamics of Earth's mantle, 
which may have start crystallizi ng at mid-depth from a magma ocean, with a surface magma 
ocean and a basal magma ocean [Labrosse et all 12007 ] . 



5.3 V + <C 1 and V ^> 1 — permeable outer boundary, impermeable inner 
boundary 

When the inner boundary is impermeable (P^ S> 1) and the outer boundary fully permeable 
(V + -C 1), the degree one mode is again found to be always the most unstable mode (Figs. [3]and 
[5]), even when 7 approaches 1. In contrast with the case where both boundary are permeable, the 
degree 1 mode now does involve deformation, and the decrease in Ra c compared to impermeable 
boundaries conditions is therefore not as drastic as when both boundary are permeable. The 
critical Rayleigh number tends toward a finite value when V + — > 0, because even with a fully 
permeable boundary, viscous dissipation always limit the development of the mode. 

This configuration could be relevant for the initiation of convection in a silicate mantle crystal- 
lizing from below from a magma ocean. The stability analysis predicts that in this configuration 
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the first unstable mode is the degree one mode shown in Figure HI However, one key point in 
this configuration is the lifetime of the magma ocean. Melting/freezing at the interface would 
play a role only if the instability growth is fast enough compared to magma ocean crystallization, 
whi ch can happen on a ky timescale in the absence of an insulating atmosphere or crystallized 
lid |Solomatovl . bood ]. 



5.4 V + > 1 and V 
boundary 



< 1 - impermeable outer boundary, permeable inner 



When the outer boundary is impermeable (T >+ 3> 1) and the inner boundary fully permeable 
(V~ <C 1), the relationship between the degree l c of the most unstable mode and 7 becomes 
non-monotonic (Figure solid blue line). l c first increases with 7, similarly to the case of 
impermeable boundaries (except l c is smaller when the inner boundary is permeable), but the 
I = 1 mode becomes again the most unstable mode when 7 exceeds ~ 0.841. Looking at the 
critical Rayleigh number as a function of I (Figure [3]), there appears to be two local minima, one 
at I = 1 and the other at a higher I, once 7 is larger than ~ 0.72. The two minima are quite close 
for all values of 7, which suggest that the I = 1 mode would be important even if it is not the 
most unstable mode. In Figure we show in blue the degree of the most unstable mode (blue 
solid line) in this configuration, as well as the degree of the local minimum at I strictly larger 
than 1 (blue dashed line). 

We show in Figured] the pattern of both the most unstable and second most unstable modes. 
The pattern of the degree one mode is found to be close to a truncated version of the pattern of 
the degree one mode of convection in a full sphere (compare with the 7 = case). 

This configuration may be relevant for the dynamics of icy satellites having an ice mantle 
overlying a global subsurface water ocean, which might be the case of several of J upite r and 



Satur n' moons, in cluding Enceladus [Nimmo and 



apnalardol. l200fil: IWaite .Tr et all . l200fll ]. Eu 



Tyler . 20081 ] . Callisto, Ganymede and Titan Spohn and Schubert , 20031 ]. Itmight also be 



ropa 

relevant for thermal convection in a silicate mantl e overlying a b a sal magma ocean, as migh t 



have been the case on Earth early in its history Labrosse et al. . 2007 ; Ulvrova et al. . 20121 ]. 



The stability analysis suggest that the length scale of convection would be significantly larger if 
melting/freezing at the interface is important. 



6 Discussion and conclusions 

The linear stability analysis presented here predicts a significant effect of semi-permeable bound- 
aries when either V~ or V + are small enough: allowing for melting/freezing at either of the 
boundary results in the emergence of larger scale convective modes. The effect is particularly 
drastic when the outer boundary is permeable, since the degree 1 mode remains the most unsta- 
ble even in the case of thin spherical shells. It seems likely that allowing for melting/freezing at 
one boundary will still result in larger scale convection at supercritical conditions, but the results 
presented here will clearly have to be supplemented by finite amplitude numerical calculations at 
supercritical conditions. In addition, the assumption of Newtonian rheology and constant viscos- 
ity limits the direct applicability of our results. The effect of variable viscosity would have to be 
investigated, in particular for application to icy moons, for which order of magnitude variations 
of viscosity across the layer may be expected. The pattern of convec t ion is also likely to depend 
on the temperature profile of the basic state McNamara and Zhong . 20051 ] . 
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At this stage, we have suggested some possible geophysical or planetological applications of 
our results, but specific studies will be needed to assess the applicability of our results in particular 
geophysical or planetological settings. In each situation, the value of V of the boundary must be 
evaluated, which necessitates some understanding of the dynamics of the melt layer in contact 
with the solid layer. 

As an example, let us discuss the case of Enceladus. Enceladus exhibits a strong hemispheri- 
cal asymme try, with the South ern hemisphere being much younger and active that the Northern 



hemisphere Porco et al.l. 1200611 . One plausible exp l anatio n for the observed asymmetry is degree 
one convec tion Grott et al. . 2007 : Steeman et al. . 20091]. Enceladus may ha ve a global subsur- 
face ocean Nimmo and Pappalardo 120061 ; Waite Jr et al. . 20091 ; Tyler . 2009 ], and it is therefore 
legitimate to consider the possible dynamical effect of melting/freezing at the inner boundary of 
the ice shell. Whether phase change at the inner boundary of the ice shell can alter s i gnific antly 



the pattern of convection depends on the value of V~: with 7 = 0.6 Schubert et all 120071 ]. the 
effect of phase change would be significant if V~ is smaller than about 10 (Figure [2]). With a 
viscosity of order 10 14 Pa s (which corresponds to the viscosity near the melting point), a radius 
R = 250 km, |A/9~| ~ 50 kg.m~ 3 and g~ ~ 0.1 m.s -2 , we find that V~ < 10 if the timescale for 
phase change is smaller than about 25 years. With given by Eq. ([7]), L = 300 kJ.kg -1 , 
c p i = 4000 K.kg _1 .K _1 , T s = 275 K, and 1 — m a ^ /m p ~ 1, thi s would require typical convective 
velocities around ~ 1 cm.s -1 in the melt layer. iTylerl 20091 ] estimates that eccentricity tides 
would have typical flow ampitude around ~ 1 mm.s -1 in a ~ 100 km thick ocean and ~ 1 cm.s -1 
in a ~ 10 km thick ocean. This would give V~ in the range 10 — 10 2 , so V~ may plausibly be 
small enough for a significant effect of melting/freezing on the pattern of convection. Including 
the effect of temperature on viscosity is likely to make the effect of melting/freezing stronger 
because the effective viscosity for relaxation of a large scale topography would be larger, possibly 
by several order of magnitude, than the high homologous temperature value of 10 14 Pa.s assumed 
here. This would yield a lower effective value of V~ , and a more permeable boundary. Whether 
or not the effect is strong enough to allow the emergence of a strong degree one convection mode 
remains an open question. The ans wer might also depend in part of the dynamical effect o f radial 
viscosity variations in the ice shell Zhong and Zuberl . l200ll ; iMcNamara and Zhong] . 120051 ] , which 
will have to be taken into account. 
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A Coefficients B[ 4 

The coefficients B\ 4 introduced in Eq. (|46l) are determined for each degree I by the boundary 
conditions at the inner and outer boundaries of the shell. 
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Using expression (pTOj) for pij , the tangential stress boundary condition (Eq. ([25]) ) gives 

B{ 1 l - 2 (l 2 -l)+B^ l l(l + 2) 
+Bi 7 - l - 3 l(l + 2) + Bi 7 - l -\l 2 - 1) 
C u( a ij7) 



(50) 



at r = 7 and 

B?(Z 2 - 1) + ^/(/ + 2) + + 2) + B{(1 2 - 1) 



C 'u( a ij) 



(51) 



at r = 1. 

The boundary condition (|25p derived from the continuity of the normal stress gives 



B{l 
1 + 



V 



I + 2 + — 7 ) + B^ 



j* — z— l 



7?- 

/ + 1 + -7 



2/ - 1 , p- 



(52) 



21(1 + 1) 



Ra 



at r = 7 and 



Bi I 1 - Z 



( / + 2 



p4 
2 

2 



3 , 1 V+\ 
- - + 1 

2 y 



21 - 1 



V 



(53) 



1 + 



a 



21(1 + 1) 



Ra 



C'u(aij) 



a 



at r = 1. 

Eqs. ([50]) . ([5T]) . ([52]) and ([53]) form a linear system of equations for l?^ 4 which is solved 
for each degree I. The B\ 4 are then used to calculate the functions Qkj in the characteristic 
equation ([17]) . 



B Translation mode 

We consider here the onset of the degree 1 mode in the limit of small Ra, V~ , and V + , but finite 
Ra/(V + +j 2 V~). With I = 1, the tangential stress free conditions ([50]) and ([5"T]) give 

^ + ^7- 5 = ^f4^ ) (54) 
Bi + Bi = Ra^l, (55) 
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while the normal stress continuity conditions (|52|) and (|53p yield 



B{V~ + 67 



+ \b{ 

T 



2 + 



a 2 7 2 



Cii(ay7) 



or 



2 + ^ 
2 



a 



(56) 
(57) 



1./ 



when "P <C 1 and V + -C 1. Using Eqs. ([54]) and ([33]) . the coefficients B 3 2 and -B3 can be 
eliminated from Eqs. (I56p and (I57f) . which give 



7 5/2 Ci 1 (a li7 ) 



2ay 

-^7^ + 3fl r f - Cil(a]J ' ) 

Noting that 

c n(«ii) = - 
Cii(ay7) 



2 ay 

2 ^|(aij7) 
2 



7raij7 



(58) 
(59) 

(60) 
(61) 



Chandrasekharl . Il96l[ ] . solving Eqs. ([35]) and ([39"]) yields 
_ J|(«ij7)/^|(«ij) -7^ Rg 

1 - 7 3/2 Ja (Qy)/Ja(aij7) ] 1 J3(ay7) 



32 



1 + 7 2 P-/T> 



1 



3 7raf - Js(aij) 

** 2 



(62) 
(63) 



It can be seen that B 2 , B^ and B\ are all ~ Ra, while i?^ ~ Ra/(V + + 7 2 7-' ). In the limit 
of small Ra, V~ ', P + , but finite Ra/{V + + 7 2 P"), we therefore have S| > To a 

good approximation, py is then given (from Eq. ([35]) ) by 



py ~ 5^r, 

and the poloidal scalar of the first unstable mode is 

P = £ Aypy (r)l?(0, 0) ~ I £ AyiJj I V Y?(0, 



(64) 



(65) 



which corresponds to a translational motion (it can be verified that a / = 1 flow with P oc r 
corresponds to a flow with uniform velocity). 

In the limit of small Ra, the characteristic equation (|47p for Z = 1 now writes 



a 



N *,k-j- S kj ~ Qkj 



0. 



(66) 
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where 



Q kj = B{ / C u (a lk r)r 5 / 2 dr. 



(67) 



Making use of recurrence relations of the Bessel functions Abramovich and Stegunl . Il965l |. we 
find that the integral in Eq. (|67[) can be written as 



C u {a lk r)r 5 / 2 dr 



ik 



Jl{aikl) 

2 



3 

7 2 



(68) 



which allows to write Qj-j as 



Qkj 



^ a lj a lk 



Js(aijj) 



J 3 {O-lj) 

Now, rewriting Eq. (|66p as 
4 



3 

7 2 



«/s(aifc7) 

2 

Js(aifc) 



3 

7 2 



(69) 



a? TVs 



-Q 



o 



(70) 



and using Sylvester's determinant theorem, we find that 



< N l,k 



Q 



Qi 



i=l 



a 2 N 3 



0, 



(71) 



from which, using Eq. (j69|) . we obtain the critical value of Ra/(V + + j 2 V ): 

-i 



Ra 



1 

- < 

4 



E 

i=l 



1 


Js(aii7) 3 1 2 
J s {an) 

2 J 


«li 


2 1 
J 2 3 (au) 

2 





(72) 
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Figure 2: Critical Rayleigh number of the I — 1 mode for : a) impermeable outer boundary, and variable 
V~ , b) impermeable inner boundary, and variable V + , c) variable V~ and V + , with V + = V~ . 
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Figure 3: Critical Rayleigh number for convection as a function of degree I, for the four end- member 
cases described in the text, for various values of 7. Note the different scales used for Ra l c . 
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(v-,r+) > 1 



(v-,r+) < 1 



v~ > 1, r+ < i 



< 1, 7>+ > 1 



Most unstable mode 2 nd most unstable mode 




Figure 4: Pattern of the first unstable mode of thermal convection in a spherical shell (streamlines), with 
aspect ratio 7 = 0, 0.2, 0.4, 0.6 and 0.8, and V~ and V + either small or large compared to 1. Impermeable 
boundaries (V >• 1) are shown by a thick line, permeable boundaries ('P ± <C 1) are shown by a thick 
dashed line. In the case V~ <C 1, V + S> 1, we also show the second most unstable mode. Only the m = 
modes are shown. 



19 




Figure 5: Degree l c of the first unstable mode as a function of the aspect ratio 7, for different configura- 
tions. The solid gray line corresponds to the case of impermeable inner and outer boundaries. The solid 
blue line corresponds to the case of impermeable outer boundary and fully permeable inner boundary. The 
dashed blue line shows the degree of the local minimum at I strictly larger than 1 in the case of imper- 
meable outer boundary and fully permeable inner boundary (see text). The solid black line corresponds 
to the cases of fully permeable inner and outer boundaries, and of fully permeable outer boundary and 
impermeable inner boundary, for which the most unstable mode is always the degree 1 mode. 
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for the translation mode in the limit of small V 



Figure 6: Critical value of the quantity (^ v++ry i V . 

and V + , as a function of the inner to ou ter radius ratio 7 (solid b lack line). The dashed black line shows 
the value of 175/2 found analytically bv lDeguen et al. submitted | for a full sphere (7 = 0). 



21 



